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1. Introduction 

A new algorithm has been developed for the Euler and Navier-Stokes equa- 
tions that uses upwind differencing based on the streamwise direction (Goorjian 
1987a, 1987b). This algorithm is time accurate and can be used in codes for 
calculating unsteady transonic flows over wings. Such codes can be used for the 
flutter analysis of wings. 

In this new algorithm, the coordinate system is locally rotated to align with 
the streamwise direction. For differencing the convective terms in the stream- 
wise direction, a new form of flux splitting is employed, in which the biasing de- 
pends cm the local Mach number. In the plane perpendicular to the stream direc- 
tion, the new flux splitting uses the condition of no flow in that local plane. (The 
formulas for the differencing in the rotated coordinate system are transformed 
to the original grid for the calculations.) By using a locally rotated coordinate 
system, the convective flux vector biasing depends on the total Mach number. 
Hence, the switching of the flux vector biasing occurs across shock waves and 
the proper domain of dependence is used in supersonic regions. For compar- 
ison, many other upwind methods switch differencing based on Mach number 
components along coordinate lines. Such criteria allow downstream influences 
in supersonic regions and switching upstream of shock waves in multidimen- 
sional flows. The formulas for the convective flux vector differencing do not 
contain any user specified parameters. Hence, the amount of numerical dissipa- 
tion is automatically determined. For viscous flows, calculations of steady flow 
over a wing using the Navier-Stokes equations were compared with experimen- 
tal data. Near the body, in a case of separated flow, the calculations showed 
improvements when compared to calculations that use central differencing with 
fourth-order dissipation terms. 

2. Flux Splitting Algorithm 

The thin-layer, Reynolds-averaged, Navier-Stokes equations in generalized 
curvilinear coordinates are given by 


d r Q + d( E + dj)F + d$G - Re 1 d$S 


( 1 ) 


2 


Goorjian 


where 


" p ' 


pU I 

pu 


puU + ZxP 

pv 

,E = J~ X 

pvU + £„p 

pw 


pwU + 6p 

. e . 


U(e + p) J 


Here Q is the vector of conservative flow variables and E is the convective flux 
vector in the £ direction. See Gooijian (1987a, 1987b) forde tails. In this paper 
the flux splitting of E will be given. The variables F and G are split in a similar 
manner. The pressure is given by 


P=( 7-1) [e-(l/2)p(u 2 + v 2 + u/ 2 )] (2) 


and U is the contravariant velocity in the £ direction. 

In solving Eq.(l) atmesh point (ijjc), the flux vector E at (i+1/2 j Jc) is split 
into E ¥ and E~, where E + is evaluated at (ij,k) and E~ is evaluated at (i+ljjc). 
Consider the case where U > 0 ; similar formulas hold for U < 0 . First the split 
flux will be given for the mass flux component E\ = J~ x pU . 

Ef = jJ-'pU [l + ( 1 - 3+)(U/q) 2 + a + (p£7/pV) 2 ] 

-s+J- x (l-l/M 2 )(p-p*)B u (q + q*)/2 (3) 

ET = jJ-'pU [l -(1 -3~HU/q) 2 -3~(pU/p*q') 2 ] 

+3-J- l ( 1 - l/M 2 )(p - p*)B u (q + q*)/2 (4) 

^ a 

The split energy flux terms and Ej are obtained by multiplying the respective 

split mass flux terms and by the local values of the total enthalpy (e+p)/p. 

The split momentum flux components in the x direction are obtained by replacing 
the mass flux terms E\ - J~ l pU in Eqs.(3) and (4) by Ei - J~ l (puU + ( x p ) » 
which is the x component of the momentum flux from E and also replacing the 
factor (q + q*)/2 by ((q + q*)/2) 2 (JQ2/p*q m ), where Q 2 - J~ x pu is the 
second component of Q. ^Similar substitutions yield the split momentum flux 
components E%, E$ and E% , E± in the y and z directions, respectively. 

The quantities p* and q * in Eqs.(3) and (4) are local sonic values of the 
density p and speed q. The quantity M is an average Mach number between 
flow quantities and their sonic values. The switches s + and s~ switch the flux 
vector biasing at sonic values, where M = 1. They equal zero for supersonic 
flow and one for subsonic flow. The quantity U is the physical component of 
the velocity corresponding to the contravariant velocity component U, where 
( C7) 2 = C/ 2 /(|V£|) 2 . The variable B u is the maximum value of p\U\/p*q*&t 
the two mesh points used to split E. For a more detailed explanation of this 
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splitting, see the papers by Goorjian (1987a, 1987b). The only changes in the 
formulas from Cartesian coordinates for curvilinear coordinates are to use con- 
travariant velocity components such as U and their corresponding physical ve- 
locity components such as U. 

3. Computed Results 

Unsteady flow 

An unsteady calculation using the Euler equations was made for flow at 
M = 0.85 over an airfoil whose thickness varies in time. Figure 1 shows the 
pressure coefficient plots for three times at which the shock wave is increasing 
in strength and moving downstream as a result of the initial thickening of the 
airfoil from zero thickness. Notice that the shock profiles at times T = 11.5 and 
T = 18.25 are sharply captured without any numerical oscillations. The shock 
profiles at those times contain one and zero mesh points, respectively. Figure 
2 shows Mach contour plots at T * 18.25. Notice that the Mach contours are 
smoothly varying in space. Also notice the tight clustering of the contour lines , 
where there is a strong transonic shock. This clustering indicates the sharpness 
of the shock capture. 

Wing C: Separated Flow 

Figure 3 shows some features of the flow field’s tip vortex. It shows traces 
of particles that were released near the tip. Figure 4 shows a comparison of 
pressure profiles at the 90% span station, which is the only station for which there 
is experimental data in the small separated flow region near the tip. At the other 
four span stations for which there is experimental data, the flow is attached and 
the two computation methods are in close agreement Notice in Fig. 4, that the 
new method, (denoted by TNS-G), produces a shock location and a suction peak 
that are in closer agreement to the experimental results than the original method, 
(denoted by TNS). The original method in the TNS code used a coefficient, DIS = 
0.2, (a relatively low value), for the fourth-order dissipation. In the region that is 
influenced by the separated flow, which is from the 90% span station outward to 
the tip, the two methods showed differences, not only in the shock wave location, 
but also along the entire upper and lower surfaces. 

Figure 5 shows Mach number contours at a height above the wing, where 
the highest Mach number is reached in the TNS calculation. Notice the higher 
contour levels reached by the TNS-G method. Also, notice that the TNS results 
show numerical oscillations in the contour levels at the downstream boundary of 
this block of the flow domain, whereas the TNS-G results do not. User specifi- 
cation of more dissipation in the TNS calculations to suppress these oscillations; 
e.g., DIS s 0.64, decreases the region of separated flow. By comparison, the 
amount of dissipation in the TNS-G results, which is automatically determined, 
suppresses the oscillations without either diminishing the separated region or 
diminishing the acceleration of the flow over the leading edge. 
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Fig 1. Pressure coefficients for the 
formation and downstream propa- 
gation of the shock wave. 


TIP VORTEX 



Fig 3. Features of the flow field’s 
tip vortex for Wing C. Plane view 
of particle traces near the dp. 



• 1.0 -.8 0 .5 1.0 1.8 


Fig 2. Mach number contours at 
T = 18.25. 
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Fig 4. Comparison of experimental 
and computed pressure coefficients 
for Wing C. 
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Figure 6 shows a comparison of boundary layer profiles at the span station, 
t; = 0 .9362 , which is near the tip, and at the chord locations x/c- 0.34 and 
x/c * 0.40, which are in the shock wave. Notice at both locations that the 
new method has a fuller velocity profile, which is another indication that near 
the body, it is automatically adding less numerical viscosity than that used with 
central differencing. 


NACA 0012 Wing: Normal Shock Flow 

Next the shock capturing features of the two methods will be compared for 
flow in which there is a shock wave separating supersonic and subsonic flow; 
i.e., a normal shock wave. On the wing surface, the pressure profiles obtained 
by the two methods are similar. However, as the flow field is examined farther 
away from the wing, the physical viscosity diminishes and the shock capturing 
features of the two calculations become determined by the differencing of the 
convective terms. 


17 * 0.93U 



Fig 6. Comparison of the boundary layer profiles for Wing C for the two com- 
putational methods. 
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For comparison, results are presented of pressure coefficients at a height of 
0.07 above the wing in units of span length. Figure 7 shows a comparison at two 
span stations. The TNS method used a value of DIS = 0.64, which is relatively 
large. Note the improvement obtained in the TNS-G results. The TNS-G results 
show no overshoots and a sharp capture of the reexpansion singularity. Also the 
sh ock is captured at the various span stations with either one or no points in the 
shock profile. For Euler calculations, these improvements are expected to occur 
at the wing surface. 


♦ TNS-G 

TNS, OIS ■ 0.64 NACA 0012 WINS 



Fig 7. Comparison of the shock capturing properties of TNS-G and TNS. Pres- 
sure coefficients in block 2 (inviscid block), at height 0.07 (in span lengths). 


References 

Goorjian, P. M. (1987a). Algorithm Developments for the Euler Equations with 
Calculations of Transonic Flows. AIAA Paper 87-0536, AIAA 25th Aerospace 
Sciences Meeting, Reno, NV, Jan. 12-15, 1987. 

Goorjian, P. M. (1987b). A New Algorithm for the Navier-Stokes Equations 
Applied to Transonic Flows over Wings. AIAA Paper 87-1121-CP, AIAA 8th 
Computational Fluid Dynamics Conference, Honolulu, Hawaii, June 9-11, 1987. 


C Cl - 

Kj Yw • • 


VJ f- 1 ; • ' “*• 

^ VI < Nl £ \ C V* L ' T ^ V*’ 


7 l 0 lO 

?i'jr 


NASA 

National Aetonautcs and 
Space Administration 


Report Documentation Page 


1. Report No. 

2. Government Accession No. 

NASA TM-101019 


4. Title and Subtitle 


A Streamwise Upwind Algorithm for the Euler and Navier-Stokes 
Equations Applied to Transonic Flows 

7. Author(s) 


Peter M. Gooijian 


9. Performing Organization Name and Address 

Ames Research Center 
Moffett Field, CA 94035 


12. Sponsoring Agency Name and Address 


National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


3. Recipient's Catalog No. 


5. Report Date 

September 1988 


6. Performing Organization Code 


8. Performing Organization Report No. 

A-88259 


10. Work Unit No. 

505-60 


11. Contract or Grant No. 


13. Type of Report and Period Covered 

Technical Memorandum 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Point of Contact: Peter M. Gooijian, Ames Research Center, MS 258-1 , Moffett Field, CA 94035 
(415) 694-6416 or FTS 464-6416 


16. Abstract 

A new algorithm has been developed for the Euler and Navier-Stokes equations that uses upwind differencing based on the 
streamwise direction (Gooijian 1987a, 1987b). This algorithm is time accurate and can be used in codes for calculating 
unsteady transonic flows over wings. Such codes can be used for the flutter analysis of wings. 

In this new algorithm, the coordinate system is locally rotated to align with the streamwise direction. For differencing 
the convective terms in the streamwise direction, a new form of flux splitting is employed, in which the biasing depends on 
the local Mach number. In the plane perpendicular to the stream direction, the new flux splitting uses the condition of no 
flow in that local plane. (The formulas for the differencing in the rotated coordinate system are transformed to the original grid 
for the calculations.) By using a locally rotated coordinate system, the convective flux vector biasing depends on the total 
Mach number. Hence, the switching of the flux vector biasing occurs across shock waves and the proper domain of depen- 
dence is used in supersonic regions. For comparison, many other upwind methods switch differencing based on Mach number 
components along coordinate lines. Such criteria allow downstream influences in supersonic regions and switching upstream 
of shock waves in multidimensional flows. The formulas for the convective flux vector differencing do not contain any user 
specified parameters. Hence, the amount of numerical dissipation is automatically determined. For viscous flows, calculations 
of steady flow over a wing using the Navier-Stokes equations were compared with experimental data. Near the body, in a case 
of separated flow, the calculations showed improvements when compared to calculations that use central differencing with 
fourth-order dissipation terms. 


17. Key Words (Suggested by Author(s)) 

Streamwise upwind algorithm 
Euler or Navier-Stokes equations 
Transonic flow applications 


19. Security Classif. (of this report) 


18. Distribution Statement 

Unlimited-Unclassified 


Subject Category: 02 


Unclassified 


NASA FORM 1626 OCT 86 


20. Security Classif. (of this page) 

21. No. of pages 

Unclassified 

7 



For sale by the National Technical Information Service, Springfield, Virginia 22161 





















